Rcpp R

[Manual] Rcpp使用手册

这是一篇关于Rcpp的基础性技术文档

Posted by Leung ZhengHua on 2017-08-30

本文总点击量次

如果你是第一次阅读这篇文章，以下资源贴可能是你感兴趣的：

Rtools下载安装
Rcpp官网
Dirk Eddelbuettel大神的神坛
Rcpp-introduction
高级R语言编程指南Rcpp章节
Rcpp Quick Reference Guide

Rcpp简介

示例一 斐波那契数列

R解决方案之一

1
2
3
4
5
fibR=function(n) {
  if(n==0) return(0)
  if(n==1) return(1)
  return(fibR(n-1)+fibR(n-2))
}

C++解决方案之一

装好Rtools最好重启一下RStudio，运行以下例子，如果编译成功，没有报错那就可以放心使用了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
library(inline)
library(Rcpp)
incltxt='
int fibonacci(const int x) {
if(x==0) return(0);
if(x==1) return(1);
return fibonacci(x-1)+fibonacci(x-2);
}
'
fibRcpp=cxxfunction(signature(xs='int'),
                          plugin = "Rcpp",
                    includes = incltxt,
                    body = '
                    int x=Rcpp::as<int>(xs);
                    return Rcpp::wrap(fibonacci(x));')

signature变量确认函数签名，定义变量进入body的接口；includes允许传递额外的指令，或者类似以下例子中的函数和类定义

R解决方案之二

前面的递归计算中，存在重复的计算，可以通过“存储”的方式更高效地计算。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
mfibR=local(
  {
    memo=c(1,1,rep(NA,1000))
    f=function(x) {
      if(x==0) return(0)
      if(x<0) return(NA)
      if(x> length(memo)) stop("你输入的值太大啦")
      if(!is.na(memo[x])) return(memo[x])
      ans<-f(x-2)+f(x-1)
      memo[x]<<-ans    #不能输错<<-，不然不能绑定局部变量
      ans
    }
  }
)
mfibR(900)

可以看到这个办法中，会在递归之前查找是否已经计算过某一个值

C++解决方案之二

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
mincltxt='
#include <algorithm>
#include <vector>
#include <stdexcept>
#include <cmath>
#include <iostream>
class Fib {
public:
Fib(unsigned int n=1000) {
memo.resize(n);  //预留n个元素
std::fill(memo.begin(),memo.end(),NAN);  //全部设为NaN
memo[0]=0.0;
memo[1]=1.0;
}
double fibonacci(int x) {
if(x<0)
return((double) NAN);
if(x>=(int) memo.size()) throw std::range_error(\"x too large for implementation\");
if(!std::isnan(memo[x])) return(memo[x]);
memo[x]=fibonacci(x-2)+fibonacci(x-1);
return(memo[x]);
}
private:
std::vector< double >memo; //存储之前的计算结果
};
'
mfibRcpp=cxxfunction(signature(xs='int'),
                     plugin = "Rcpp",
                     includes = mincltxt,
                     body='
                     int x= Rcpp::as<int>(xs);
Fib f;
return Rcpp::wrap(  f.fibonacci(x-1) );
                     ')
mfibRcpp(600)

R解决方案之三

1
2
3
4
5
6
7
8
9
10
11
12
fibRiter=function(n) {
  first=0
  second=1
  third=0
  for (i in seq_len(n)) {
    third=first+second
    first=second
    second=third
  }
  return(first)
}
fibRiter(1000)

由于迭代算法既不需要存储状态，也不需要递归，所以其可以在存储方法的基础上进一步提高。

C++解决方案之三

1
2
3
4
5
6
7
8
9
10
11
12
13
14
fibRcppIter=cxxfunction(signature(xs="int"),
                        plugin = 'Rcpp',
                        body='
                        int n=Rcpp::as<int>(xs);
double first=0;
double second=1;
double third=0;
for(int i=0;i<n;i++) {
third=first+second;
first=second;
second=third;
}
return Rcpp::wrap(first);
                        ')

这是目前最快的版本。

数据结构

IntegerVector 类

返回完美数

完美数是指等于其因子之和的正整数，本例的函数签名为空，选择“Rcpp”作为plugin会引导cxxfunction()寻找Rcpp中的头文件和库：

1
2
3
4
5
6
7
8
9
10
src='
Rcpp::IntegerVector epn(4);
epn[0]=6;
epn[1]=14;
epn[2]=496;
epn[3]=8182;
return epn;    //这里没有显式声明从C++向量转换为R向量
'
fun=cxxfunction(signature(),src,plugin = 'Rcpp')
fun()
  • 由于隐式调用wrap，返回向量不需要额外的代码，直接使用return epn;
  • 创建一个新向量和定义其大小可以调用Rcpp::

使用R向量输入

1
2
3
4
5
6
7
8
9
10
src='
Rcpp::IntegerVector vec(vx);
int prod=1;
for (int i=0;i<vec.size();i++) {
prod*=vec[i];
}
return Rcpp::wrap(prod);  //这里的返回值用warp包裹起来了
'
fun=cxxfunction(signature(vx='integer'),src,plugin = 'Rcpp')
fun(1:10)

NumericVector类

使用两个输入

1
2
3
4
5
6
7
8
9
10
11
12
src='
Rcpp::NumericVector vec(vx);
double p=Rcpp::as<double>(dd);
double sum=0.0;
  for(int i=0;i<vec.size();i++) {
    sum+=pow(vec[i],p);
  }
return Rcpp::wrap(sum);
'
fun=cxxfunction(signature(vx='numeric',dd='numeric'),src,plugin = 'Rcpp')
fun(1:4,2.2)
FEATURED TAGS
Rcpp R
FRIENDS